Chaotic stirring by a mesoscale surface-ocean flow 
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The horizontal stirring properties of the flow in a region of the East-Australian Current are 
calculated. A surface velocity field derived from remotely sensed data, using the Maximum Cross 
Correlation method, is integrated to derive the distribution of the finite-time Lyapunov exponents. 
For the region studied (between latitudes 36°S and 41°S and longitudes 150°E and 156°E) the mean 
Lyapunov exponent during 1997 is estimated to be Aoo = 4x 10~ 7 s _1 . This is in close agreement 
with the few other measurements of stirring rates in the surface ocean which are available. Recent 
theoretical results on the multifractal spectra of advected reactive tracers are applied to an analysis 
of a sea-surface temperature image of the study region. The spatial pattern seen in the image 
compares well with the pattern seen in an advected tracer with a first-order response to changes in 
surface forcing. The response timescale is estimated to be 20 days. 
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Satellite imagery of the ocean surface shows that 
sea-surface temperature and chlorophyll have 
complex distributions. The filaments and whorls 
characteristic of tracers in stirred fluids are of- 
ten clearly evident. Despite the importance of 
stirring to a range of problems, such as plank- 
ton ecology, larval transport or predicting the 
fate of pollutants, it has not been well described. 
Here we analyse the stirring properties of a ve- 
locity dataset obtained from the interpretation of 
sequential satellite images. The finite-time Lya- 
punov exponents of the flow are determined. Us- 
ing these results, we show that the patterns seen 
in a sea-surface temperature image are consistent 
with those formed from the advection of a reac- 
tive tracer. This suggests that an understanding 
of the distribution of tracers such as tempera- 
ture and chlorophyll in the surface ocean must 
include a representation of both the chaotic ad- 
vection by the flow and the dynamics of the trac- 
ers themselves. If the advection can be charac- 
terised, as we have been able to do here, then it 
may be possible to use the information contained 
in sea-surface imagery to infer the time-scales of 
the tracer dynamics. 
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INTRODUCTION 



One of the simplest measures of the distribution of a 
tracer is its power-spectrum, which quantifies the vari- 
ability in the distribution at each spatial scale. Between 
inverse wavenumbers of 1 and 100 km the spectra of sea- 
surface temperature and chlorophyll are usually found to 
have a power-law form, F{k) ~ . The exponent, (3, 
is typically in the range 1.5 < f3 < 2.5 @, |, |, |, § §. 
Where both sea-surface temperature and chlorophyll are 
measured at the same time their spectra often have sim- 
ilar slopes. Ocean currents are turbulent within the 
mesoscale wavenumber range, 1 km < fc _1 < 100 km, 
and their energy falls off rapidly at higher wavenum- 
ber J?j]. The flow within this range is approximately 
two dimensional, and it is appropriate to consider the 
advection of a passive tracer in the surface ocean as a 
stirring process ]8|, M. Stirring acts to stretch and fold 
patches of tracer, distorting them into long tendrils and 
filaments. In spectral terms, variance that is input at 
low wavenumber is transferred by stirring towards higher 
wave-numbers, where it is dissipated by diffusion. The 
spectrum of a stirred tracer, which has a source at low 
wavenumber but is otherwise conserved, is expected to 
be power-law with an exponent /? = 1 |lCj] . Batchelor 
law spectra with this form are sometimes seen in labo- 
ratory experiments of stirring by two-dimensional flows 
pH , but (3 = 1 is outside the range usually measured 
for sea-surface temperature and chlorophyll, as has been 
noted by several authors jj], ||. It is clear, however, that 
neither sea-surface temperature or chlorophyll are con- 
served quantities. For a non-conserved tracer which sat- 
isfies a first-order equation, any spectral exponent /3 > 1 
can be obtained, depending on the relation between the 
reaction-rate and the typical stirring-rate of the flow [|l2j . 
This suggests that the patterns seen in the surface ocean 
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may reflect both the stirring by the flow and the tracers' 
non-conservation. Recently, results have been derived 
which give the multifractal structure of an advected re- 
active tracer in terms of the distribution of finite-time 
Lyapunov exponents, which characterise the stirring by 
the flow [[[3) . The theoretical results were derived in the 
context of Lagrangian chaotic flows. Applying the theory 
to the surface ocean is hampered by the lack of data on 
stirring by ocean currents. The few measurements which 
have been made depend on following the stretching of a 
tracer patch for an extended period of time |TJ, |TJ, [l6| . 
This technique gives little information on the distribu- 
tion of stirring rates, each observation of a tracer patch 
typically returning only a single estimate of the stirring. 
Given the lack of direct measurents, an indirect approach 
must be taken. In this paper, surface velocity data, in- 
ferred from satellite observations, are used to calculate 
stirring rates. The analysis is focussed on the East Aus- 
tralian region, where there is strong mesoscale eddy ac- 
tivity and where a suitable dataset is available. The de- 
rived distribution of stirring rates, more formally finite- 
time Lyapunov exponents, is applied to an analysis of the 
distribution of sea-surface temperature. 
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SURFACE VELOCITY DATA 
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It is possible to estimate sea-surface velocity from a 
comparison between sequential satellite images of sea- 
surface temperature jl7|, [l8|, 19 . This Maximum Cross- 
Correlation (MCC) method has recently been applied to 
7 years of satellite data from 1993 to 1999 pC] , cover- 
ing the region where the East Australian Current (EAC) 
separates from the coast and heads across the Tasman 
Sea. The MCC data has a high spatial and temporal 
resolution, but because no sea-surface temperature im- 
agery can be obtained through cloud, there are extensive 
data dropouts. This problem has previously limited the 
application of the MCC method to sequences of a few 
exceptional cloud free images. In the dataset used here, 
the MCC velocities have been blended with lower resolu- 
tion data derived from an analysis of satellite-altimeter 
measured sea-surface height [^l] . The resulting optimally 
interpolated (01) surface velocity field has a uniform spa- 
tial and temporal resolution of 100 km and 10 days, suffi- 
cient to resolve mesoscale features. While both the satel- 
lite altimeter and temperature data are globally avail- 
able, they have only been combined to derive surface ve- 
locities in the region that is studied here. In this paper 
the stirring properties of the OI flow are characterised. 

The EAC current is an intense western boundary cur- 
rent which flows south, close to the Australian coastline, 
until it reaches the latitude of Sydney. Here the current 
turns, forming large loops that pinch off to form eddies. 
The eddies continue moving south while the main current 
turns to flow east. The separation region is seen as a peak 
in the eddy kinetic energy of the flow (Fig. |l|) . Analysis 
is restricted to a smaller region to the south of the main 



FIG. 1: The study region in the south-west Tasman Sea, 
showing the region over which the satellite derived velocity 
data is available. The contours show the average eddy kinetic 
energy (half the velocity variance, cm 2 s -2 ) from the 1997 
data. The small box marks the area over which stirring rates 
are derived. 



separation zone, between latitudes 36°S and 41°S, and 
between longitudes 150°E and 156°E. The mean currents 
in this region are relatively weak, so trajectories remain 
in the area where the OI velocities are defined for long 
enough to allow the stirring to be determined. 

In the following two sections, results on stirring in two- 
dimensional flows are briefly reviewed. 



III. FINITE-TIME LYAPUNOV EXPONENTS 

In a two-dimensional divergence-free flow, with veloc- 
ity v = (it, v), the change in the separation between the 
trajectories of two infinitesimally separated points x(t) 
and x(t) + Sx(t) is 



Sx = S(t)Sx, 
where the matrix S(t) is the Jacobian, 



S(t) 



d x u d y u 
d x v d v v 



(1) 



(2) 



and is evaluated along the trajectory x(t) [Q. The Ja- 
cobian may be resolved into a symmetric and antisym- 
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metric part. The antisymmetric part is a pure rotation. 
The symmetric part is the strain tensor, 



E(t) 



d x u (d y u + d x v)/2 
(d y u + d x v)/2 -d x u 



(3) 



which acts to stretch a small patch of tracer without 
changing its area. The strain tensor has the eigenvalues 
7± = ±yjd x u 2 + (d y u + d x v) 2 /4: 7 with the correspond- 
ing eigenvectors 'y±. In a pure-strain flow with velocity 
v = ("yx, —jy) the strain tensor is diagonal, and the solu- 
tion to Eq. ((1) is Sx = (fe(O)e 7 *, 5y(0)e" 7 *). There is an 
exponential growth in the component of the separation 
which is aligned with the positive eigenvector, and an 
exponential decay of the other component. In this flow 
field a small initially circular patch of tracer is stretched 
into an ellipse, with the rate of growth of the major axis 
being given by the strain-rate, 7. 

In a general time-dependent flow infinitesimal separa- 
tions between particles will still grow exponentially with 
time, but the growth is no longer a simple function of 
the strain |^| . It may be characterised by the finite-time 
Lyapunov exponents, defined by 



(4) 



where the orientation of the initial separation, 5x (0), is 
chosen so that A is maximal. At large times, A may 
be approximated by using Eqs. (0) and (@) wi * h 

a ran- 
domly chosen initial orientation of 5 the stretching 
by the flow aligns most initial vectors with the direction 
of maximal stretching. We refer to the Lyapunov expo- 
nents approximated in this way as A. To calculate the 
Lyapunov exponents at small times more care must be 
taken. An infinitesimal circle is deformed by the flow into 
an ellipse. The integrated deformation along a trajectory, 
M(t), may be obtained by solving 



M = SM, 



(5) 



where M(Q) is the identity matrix. The semi-major and 
semi-minor axes of the ellipse are the eigenvectors of the 
matrix M T M, with their squared lengths being the cor- 
responding eigenvalues (denoted m + and to_, with m+ 
referring to the larger value). The finite-time Lyapunov 
exponent along a trajectory is then given by 



X(x(t),t) = — log(m+), 



(6) 



In a closed ergodic flow the Lyapunov exponents con- 
verge with time towards a single value, lim^oo A(t) = 
Aoo. If the flow is chaotic then > 0. At times 
t >> 1/Aqo the averaged finite-time Lyapunov exponent 
has the form 



(7) 



(A) - A/t + B/Vi+ A 



where A and B are two constants p4| . This relation 
has previously been used to estimate Aoo where the in- 
tegration time was too short to allow convergence of the 



finite-time exponents to be achieved [g5| . The long-time 
behaviour of the standard deviation of the finite time 
Lyapunov exponents, a\, is given by 

a x ~ y/A/t, (8) 

where A is a third constant ]24j| . The probability distri- 
bution of finite-time Lyapunov exponents has the time- 
asymptotic form 



P(A(t),i)^^ 2 e - G ^, 

where G(A) > and G(A oc ) = G'iX^) = 0@ 
P(X,t) is gaussian then G(A) is a parabola, 



G(A) 



(A~Aoo) ; 
2A 



(9) 
2§. If 

(10) 



IV. 



THE DISTRIBUTION OF ADVECTED 
TRACERS 



As a simple approximation, sea-surface temperature 
may be regarded as dynamically passive, with horizontal 
temperature gradients not affecting the flow. This view is 
supported by observations which show that, away from 
major fronts, horizontal temperature variations within 
the surface layer are matched by salinity variations in 
such a wa y t hat there are only weak horizontal density 
gradients [g7| £29) . As a parcel of surface water is ad- 
vected by the flow it exchanges heat with the atmosphere 
and with the underlying water. An initial exploration of 
the dynamics of sea-surface temperature may be made 
by assuming that the changes in temperature are first- 
order. Consider a tracer C(x(t),t) which satisfies the 
Lagrangian equation 



dC 
~dt 



a(Go-G), 



(11) 



where a is a relaxation rate. The source, Co(x,t), is 
assumed to only vary over large scales, i.e. there is a 
wavenumber fco such that the fourier power spectrum, 
|Co(fc)| 2 , is zero for k > k . If the tracer is advected by 
a flow which has X(t) — Xqo and which has a horizontal 
diffusivity D, then for wavenumbers between fco and the 
diffusive cut-off, fc n = y/ Xqq/D, the power spectrum has 
a power-law form Jl2|, 



\C{kf\ 

The power-law exponent is 



(12) 



(13) 



For tracers with a rapid relaxation time there is less 
transfer of variance toward higher wavenumber and the 
power-spectrum is steep, whereas for tracers with a slow 
relaxation rate the spectrum approaches that of a forced 
conserved tracer with Batchelor law scaling, (3—1 JlOfl . 
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For more general flows with a distribution of finite- 
time Lyapunov exponents (Eq. ^J) the tracer distribution 
is multifractal, with the scaling exponents being related 
to the function G(A). Multifractal analysis has previ- 
ously been carried out on small-scale temperature and 
phytoplankton data (e.g. |3(], g^V. Here theoretical re- 
sults obtained by Neufeld et al. p3| are presented. The 
qth order structure function is defined as 



S q (Sr) 



(\C(x + 6x,t)~C(x,t)\ q } 



(14) 



where Sr = \Sx\ and the angle brackets indicate averaging 
over spatial locations and separations Sx. For a scale- free 
distribution, the structure functions have a power-law 
form in the limit Sr — > 0, 



S q (6r) ~ St**. 



(15) 



For an advected first-order tracer (Eq. y_l|) the scaling 
exponents are given by 



( q = min <^ q 



<!<> ■ A, j 



A 



(16) 



where the minimum is taken over all values of A > 0. 
The power-spectral exponent is related to the structure 
functions by the relation 



= Ca + 1. 



(17) 



In cases where the function G(A) may be approximated 
by a parabola (Eq. |l0|) the scaling exponents have the 
form 



A, 



C, = min iq, ^(Vl + 2qaA/\%, - 1) j . (18) 

In the limit A — > 0, P(X) becomes a delta function, S(X — 
Aoo), and the tracer distribution is monofractal with the 
scaling exponents 



C q = qa/ Aoo (for a < Aoo). 



(19) 



When q = 2 the equation for the power-spectral expo- 
nent, Eq. (|l3|), is recovered using Eq. (|l~7"| ). The scaling 
exponents have the same form (Eq. |19[) in the q — > 
limit. This may be used to estimate a/Xoo even if the 
function G(A) is unknown. 



V. FINITE-TIME LYAPUNOV EXPONENTS OF 
THE OI DATA 

The finite-time Lyapunov exponents are calculated for 
the study region using the 1997 OI velocity data. A 4th 
order Runge-Kutta integration with a daily time-step is 
used to calculate the trajectories of a regular grid of 
points. The OI velocities are objectively mapped onto 
a regular grid with a 20 km by 20 km by 5 day spacing, 
and the data are linearly interpolated in time and space 



to estimate the velocity along the trajectories. The Jaco- 
bian is integrated following Eq. (^]) to obtain the matrix 
M T M, and the Lyapunov exponents are then calculated 
from the eigenvalues m + (Eq. ||). There are two compo- 
nents to the OI velocity field, the mean flow and the ve- 
locity anomaly. The optimal interpolation ensures that 
the anomaly is divergence free, but the mean flow is de- 
rived directly from the MCC data and does not have this 
constraint. While the full velocity field is used to derive 
the trajectories, only the traceless component of the Ja- 
cobian is integrated to obtain M . As an example of the 
results, the flow is shown on the 22 November 1997 in 
Fig. |^, along with the Lyapunov exponents for three dif- 
ferent times t. In this figure a regular grid of points, with 
a spacing of 2.2 km, was integrated backwards in time. 
Trajectories that left the region within which the OI ve- 
locities were defined are shown in black. The Jacobian 
was then integrated forwards along the trajectories to ob- 
tain the Lyapunov exponents at a regular grid of points. 
The streamfunction on November 22 is dominated by an 
anti-cyclonic eddy with a diameter of approximately 200 
km, centred on 38°S 152°E. The ten day Lyapunov expo- 
nents show that the eddy centre is an area of relatively 
low stretching, with the highest stretching rates being 
in arms around the eddy. For larger times t the Lya- 
punov exponent develops a filamentary structure, with 
the width of the filaments narrowing with time. 

The probabability distribution P(X,t) is shown in 
Fig. H for various times t. This was calculated from sets 
of forward trajectories initialised on a 61 x 51 regular 
grid. Each set of trajectories began 10 days apart, be- 
tween 1 January and 12 October 1997. As expected, the 
distribution of the exponents narrows with time as each 
trajectory experiences a range of flow conditions. The in- 
tegration is not continued beyond 80 days as by that time 
over a third of the trajectories have left the region where 
the OI velocities are available. The strain of the OI flow 
in the study region has the mean and standard deviation 
(7+) = 2.5 xl0~ 6 s- 1 andcr 7+ = 1.6 x 10~ 6 s -1 . The ini- 
tial stretching is at the rate given by the strain, but the 
average Lyapunov exponent decreases with time (Fig. |4|). 
The mean Lyapunov exponent is still decreasing at the 
end of the 80 day integration so the form of G(A) cannot 
be determined. The Lyapunov exponent Aoo can only be 
estimated by fitting the relation for the expected time- 
dependence, Eq. (Q) to the data. The least-squares fitted 
curve to the mean Lyapunov exponent between 15 and 
80 days is (A) - 0.371/4 + 4.17 x IQ^/^ft + 3.84 x 10~ 7 , 
where the time, t, is in seconds. This gives the estimate 
Aoc = 4 x 10" 7 s" 1 (or Aoo = 0.035 day^ 1 ). A similar fit 
of Eq. (||) to the standard deviation of the Lyapunov ex- 
ponents gives A = 5 x 10 -7 s _1 . We have not attempted 
to quantify the uncertainties of these estimates, but, be- 
cause of the slow convergence, the uncertainties will be 
relatively large. It could be expected that P(X, t) would 
become gaussian when a\ « A^, or when t >> A/A^.. 
For the OI flow this will occur when t >> 40 days, so 
it is clear that the integration is too short to allow the 
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FIG. 2: The streamfunction and the finite-time Lyapunov exponents within the study region on 22 November 1997. (a) the 
streamfunction, ip (xlO 4 m 2 s" 1 ) (b) X(t), t = 10 days (c) \(t), t = 20 days (d) \(t), t = 40 days. 
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FIG. 3: The probability distribution of the finite-time Lya- 
punov exponents. As expected the distribution narrows with 
time, but even after 80 days the distribution has not con- 
verged. 



asymptotic form of P(A, t) to emerge. 

The mean value of the Lyapunov exponents approx- 
imated from the stretching of an initially randomly 
aligned vector, A, is shown in Fig. |] for comparison. 
Because of the random alignment of the initial vector 
relative to the straining of the flow, the initial stretch- 
ing is zero. It grows rapidly, and then slowly con- 
verges toward the actual value, (A). A least squares 
fit of Eq. (0) to the data between 15 and 80 days gives 
(A) - 0.0428/i + 1.59 x lO^/V* + 4.33 x 1CT 7 . This is 
consistent with the estimate of A~, — Ax 10~ 7 s~ 1 . 



VI. SEA-SURFACE TEMPERATURE 

In this section the formalism developed by Neufeld et 
al. p3|] and presented in section IV is applied to a satel- 



lite image of sea-surface temperature (SST) to test the 
applicability of the theory in this context. In Fig. ^| an 
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FIG. 4: The variation with time of the mean (circles) and 
standard deviation (squares) of the distribution of finite-time 
Lyapunov exponents. The open symbols are from the distri- 
bution of the maximum Lyapunov exponents, A. The solid 
symbols are from the approximation to the Lyapunov expo- 
nents, A, calculated from the stretching of randomly aligned 
initial vectors using Eq. (Q). The solid lines show the least- 
squares fit to the data between 15 and 80 days, using Eqs. ([?]) 
and (^). The horizontal dashed line marks the asymptotic 
value of the Lyapunov exponent, Aoo = 3.84 x 1CP 7 s _1 . 



SST image from 22 November 1997 is compared with 
images produced by the advection of a first-order tracer 
by the OI flow. The OI streamfunction for this day is 
shown in Fig. |2|, with the dominant feature being an anti- 
cyclonic eddy. The equation for the tracer (Eq. [ll]) is 
integrated along the trajectories, which were calculated 
backwards from a regular grid with a 2.2 km spacing. 
The integration is only for 40 days as for longer times 
too many trajectories leave the area in which the OI ve- 
locities are defined. The function Cq is taken to be a 
linear gradient given by the latitude of the trajectory. 
The tracer field along each trajectory is initialised at its 
starting latitude. With a relaxation time of 1/q = 20 
days (Fig. ||(b)) the broad structures of the modelled 
SST compare well with what is seen in the satellite im- 
age taken on that day (Fig. |](a)). The anti-cyclonic eddy 
appears as an area of low gradient, with the northern 
boundary being marked by a line of strong gradient at 
37° S. A tongue of warmer northern water comes down 
the western side of the eddy, and there is a meander- 
ing front between latitudes 39°S and 40°S. Given the 
extreme simplicity of the sea-surface temperature model 
the agreement between these images is striking. There 
are, however, many smaller features which are not re- 
solved in the OI velocity field and so the contours of the 
modelled SST are smoother than the contours of the ac- 
tual SST image. If the relaxation rate is set to a — 0, and 
the latitude is advected as a conserved tracer, then the 
same broad patterns are seen (Fig. ||(c)), but the tracer 
structure becomes more striated than the observed sea- 
surface temperature. The highest tracer gradients are in 
the regions of greatest stretching, and the filaments of 




Longitude (°E) 

FIG. 5: A satellite sea-surface temperature image from 22 
November 1997 (a) is compared with SST modelled using 
the OI velocity field and a first-order equation for a reac- 
tive tracer, Eq. ([n]). The two modelled SST images were 
generated with the relaxation rates (b) a = 0.05 day -1 and 
(c) q = 0. In (a) black marks either cloud or invalid data, 
and in (b) and (c) black marks trajectories which went out- 
side the area in which the OI velocities are defined. The SST 
image has 1 km pixels, and the modelled data has a 2.2 km 
resolution. 
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FIG. 6: Power spectra of the actual and modelled SST data 
shown in Fig. ^. The spectra are the mean of the 1-D power 
spectra taken along 250 km long sections of constant latititude 
which contain no invalid data. The normalisation is arbitray, 
having been chosen to aid comparison of the spectra. The 
spectrum of the actual SST data is shown by the line with 
circles. The other lines correspond to the modelled SST, with 
the following values of a (day -1 ), from the steepest to the 
flattest: 0.2, 0.1, 0.05, 0.025 and 0. 



the t — 40 day Lyapunov exponents (Fig. ||) are aligned 
with the contours of the tracer (Fig. |](c)). 

The spatial patterns in these images may be quanti- 
tatively compared through their power-spectra, Fig. [(]. 
The spectra were taken from the 250 km long constant- 
latitude sections through the image which only contained 
valid data. The spectra shown are the averaged cyclic 
power-spectra, with the sections being detrended and 
multiplied by a Hann filter before calculating the peri- 
odogram. The averaged power-spectra have been given 
an arbitrary normalisation to allow them to be readily 
compared. Between inverse wavenumbcrs of 10 and 100 
km the power spectra have a power-law form, with there 
being a progressive steepening of the spectra as a in- 
creases. This has previously been seen in simulations of 
a first-order tracer in a model flow intended to represent 
mesoscale turbulence p2[ , and is as expected from the- 
ory (Eq. [F^) . The power-law exponents of the spectra in 
Fig. | are = 2.44 ± 0.04 (actual SST); = 2.12 ± 0.03 
(conserved tracer); = 2.28 ± 0.03 (1/a = 40 days); 
= 2.46 ± 0.04 (1/a = 20 days); = 2.77 ± 0.09 (1/a = 
10 days); and = 3.47±0.14 (1/a = 5 days). The spec- 
tral exponent of the actual SST data is not significantly 
different from the exponent of the modelled data with 
1/a = 20 days, suggesting that the SST is responding to 
changes in forcing with a timescale of around 20 days. 

The spectrum of the conserved tracer is flatter than 




FIG. 7: Structure functions, S q (Sr), calculated using 
Eq. @ from the November 22 SST image (Fig. |(a)). The 
scaling exponents, £((/), are calculated from a least-squares 
fit to the structure functions within the range 5 to 30 km, 
shown by vertical dotted lines. The structure functions are 
approximately power-law over this range, but roll off toward 
larger separations, Sr. 



the tracer with 1/a = 40 days but, because of the short 
integration time, the spectrum is still steeper than the 
k^ 1 form expected for a forced conserved tracer. The 
40 day period over which the trajectories are defined 
is too short to allow a full transfer of variance from 
low to high wavenumbers. The timescale expected for 
the transfer between two lengthscales L max and L m j n is 
T L = log(L raax /L min )/A oc . For the OI flow 1/Aoo = 29 
days, so with £ max /L m j n = 10 the timescale for variance 
transfer through the mesoscale is Tr, = 70 days. An inte- 
gration of at least this length would be needed to allow a 
spectrum with a slope close to to be established. The 
40 day trajectory length is too short to allow the relation 
between spectral-exponent and the relaxation-rate, given 
in Eq. ([l3|), to be tested. 

The multifractal scaling of a tracer may be calculated 
using the equation for the structure function (Eq. fl4"| ). 
The structure functions of the SST image shown in 
Fig. |(a) and of simulated SST were calculated for Sr < 
50 km, and with q between 0.2 and 3 (Fig. ^). The mod- 
elled data were calculated as for Fig. |5|(b, c), but with an 
80 day integration time and with a range of relaxation 
rates a. The longer trajectory lengths meant that 25% 
of the trajectories went outside the area where the ve- 
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FIG. 8: Multifractal scaling exponents of sea-surface temper- 
ature. The solid stars (*) are the exponents calculated from 
the SST data shown in Fig. |H](a). The solid line is a least- 
squares fit of Eq. @ to the exponents. The other symbols 
and the dotted lines mark the scaling exponents calculated 
from the modelled data and the associated best-fit curve, with 
the following values of a (day^ 1 ): 0.2 (circle), 0.1 (square), 
0.05 (diamond), 0.025 (triangle), (plus). 



locities were defined. Despite the presence of the invalid 
data, it was still possible to calculate the structure func- 
tions as the averaging is only taken over pairs of valid 
points. The data were detrended before the calculation 
was made. The resulting scaling exponents, obtained 
from a least squares regression of \og(S q (Sr)) against 
log(Jr), are shown in Fig. |L The structure functions 
of the SST image lie close to the curve derived from the 
simulated SST with 1/a = 20 days, in agreement with 
the estimate of a obtained by comparing power-spectra. 
A least-squares fit of the relation given by Eq. (|l^) closely 
follows the SST data, with the best fit dimensionless pa- 
rameters being a/Aoo — 0.90 and A/Aoo = 1.12. With 
the value estimated for the OI flow of Aoo = 4 x 10 -7 s _1 
this results in a — 0.031 day~\ or 1/a = 32 days, and 
A = 4.5 x 10~ 7 s _1 . There is a long chain of assumptions 
required to allow this estimate to be made, but it doesn't 
require an explicit modelling of the flow or, indeed, any 
knowledge of the flow beyond Aoo. Although it is perhaps 
coincidental, it is pleasing that the value of A is close to 
that found from an analysis of the finite-time Lyapunov 
exponents. This suggests that the spatial structure seen 
in the satellite image is consistent with the representa- 
tion of SST as an advected first-order tracer. The other 
dotted lines in Fig. ^ result from fitting the expression for 
((q) to estimates made from the simulated SST. The val- 
ues of 1 1 a (days) predicted from the best-fit curves are 21 
(5); 24 (10); 29 (20) and 50 (40), where the actual value 
is shown in brackets. The curve fitting over-estimates 
the relaxation time-scale, particularly when the relax- 
ation is rapid. This is likely to be because the derivation 
of the expression for the fitted curve (Eq. |l8|) assumes 
that P(A, t) is gaussian. This is a poor approximation at 



small times. 

Is a 20 day response time for sea-surface temperature 
reasonable? A one-dimensional turbulence closure model 
was recently used to simulate seasonal variability of SST 
in the eastern Tasman Sea, at a similar latitude to the 
data analysed here, but close to New Zealand |||. The 
model used heat and wind forcing derived from meteo- 
rological observations. It was found that a perturbation 
of the atmospheric forcing, applied so that the perturba- 
tion was removed at the end of November, had an effect 
which decayed with a time-scale of one to two months. 
The response was more complex than a simple first or- 
der relation, but nevertheless the response-time inferred 
from the water column model was of a similar order to 
that found here. It would be interesting to combine the 
Lagrangian advection with a water-column model that 
included realistic heat and wind forcing. 

A similar set of assumptions about the dynamics of 
sea-surface temperature has been used in mesoscale sim- 
ulations of SST variability by Klein and Hua (3^]. The 
focus of these simulations was on the response of SST 
to episodic bursts of wind. When the mixed-layer is 
deepened by a wind event the surface temperature is im- 
printed with a signal from the sub-surface ocean. The 
mesoscale dynamics below the thermocline may be repre- 
sented by quasi-geostrophic turbulence. In the ocean in- 
terior, temperature is not passive and is dynamically con- 
strained, with a predicted spectral exponent of (3 = 3 in 
the sub-mesoscale range. So, when the mixed-layer deep- 
ens the sea-surface temperature spectra steepen towards 
this value. During subsequent stirring of the mixed-layer 
it was found that temperature variance was transferred 
towards higher wave-number, resulting in the spectral 
exponent becoming smaller with time. From this point 
of view, the processes determing sea-surface temperature 
spectra are a steepening through episodic wind-mixing 
and a flattening through horizontal stirring. Following 
this argument, sea-surface temperature spectra with ex- 
ponents less than one should be possible during periods 
of time when the mixed-layer is not deepening for many 
months, such as over summer. Such flat spectra are not 
usually observed. An analysis of the spatial pattern in 
a time-series of SST imagery would allow the relative 
importance of continuous relaxation and episodic wind- 
mixing to be resolved. 



VII. DISCUSSION 

Estimates or descriptions of stirring in the surface 
ocean are rare. In this paper, a spatially and tempo- 
rally complete velocity dataset was used to calculate the 
Lyapunov exponents of the flow in a region of the East- 
Australian current. The mean stirring rate was found to 
be Aoo = 4 x 10~ 7 s _1 , corresponding to a timescale of 
1/Aoo = 29 days. The velocity field used is derived us- 
ing relatively indirect techniques, and only resolves the 
largest mesoscale features. The question then arises, how 
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reliable is the estimate of the stirring? The observa- 
tion that tracers advected in ocean flows are filamental 
suggests that there is a separation between the diffusive 
length scale and the scale which is controlling the stirring. 
If, indeed, the larger mesoscale features are dominating 
the shear then the analysis will be appropriate. This 
will not be the case in shallow waters, where the eddy 
length scale is too small to be resolved using the remote 
methods that are relied on here. There have been beau- 
tiful observations made of sub-mesoscale stirring, seen in 
photographs taken from space of sun-glitter on the sea- 
surface [^5[ ||(|. These show small eddies, with diameters 
of only tens of kilometers, organising surface flotsam into 
filamental slicks. It appears from these photographs that 
sub-mesoscale eddies may at times dominate the stirring 
of tracers. In this case, the methods used here would 
underestimate stirring rates. 

From the distribution of the Lyapunov exponents it is 
clear there is a very wide range in the stretching rates 
which could be experienced by a patch of tracer in the 
ocean. Moreover, it would be difficult to predict the 
stretching of a patch from an Eulerian snapshot of the 
flow. While the analysis used here is based on a veloc- 
ity field which could, in principle, be obtained in near 
real time, predicting the dispersal of either a pollutant 
or a deliberately released tracer is unlikely to be possi- 
ble. Small errors in the velocity will lead to a rapid di- 
vergence of trajectories. Despite the spread of the prob- 
ability distribution of Lyapunov exponents it is intrigu- 
ing to note that the only previous estimate of stirring in 
the surface ocean, obtained from a satellite image of a 
deliberately induced phytoplankton bloom 16 , found a 
value for the stretching of A(42 days) = 6 x 1CP 7 s _1 . 
Other estimates of mesoscale stirring from a numerical 
model (5.8 x 1CP 7 s _1 |Bj) and from a deep tracer re- 
lease (3±0.5 x 1(T 7 s _1 pf H)) are all of a similar order. 
The agreement of the values obtained, using very differ- 
ent methods, suggests that the estimate for Aoo found by 
analysing the 01 velocity field is a reasonable one, with 
the advantage over the tracer-releases that the required 
data are readily available. 

There has been much work recently on the stirring of 
reactive tracers, particularly phytoplankton, in the sur- 



face ocean || H H U El Hi H @ ©• This re- 
search has generally relied on simple models of the ocean 
flow. The study of Lagrangian chaotic flows, in partic- 
ular, has allowed the problem to be broken into its sim- 
plest components and understanding of the patterns seen 
in advected tracers has been advanced considerably. The 
application of this theoretical work has been somewhat 
hampered by the lack of appropriate oceanic data. The 
analysis presented here is intended to help alleviate that 
problem, but is presented as a simple case study only. 
There are many further directions in which it may be 
taken. It would be interesting to test the application 
of the theoretical relationships in a model of ocean tur- 
bulence which was not disadvantaged by the short inte- 
gration times that we were restricted to here. As far as 
sea-surface temperature is concerned, it would be inter- 
esting to carry out a systematic analysis of the whole 
East- Australian optimally- interpolated dataset. There 
may be seasonal variations in the estimated value of the 
response time, a, and if these reflect real ocean processes 
then they will contain information on mixed-layer depth. 
Since mixed-layer depth is a crucial factor controlling pri- 
mary production in the ocean, a method for deriving it 
from remotely sensed data would have great value. As 
more satellite data becomes available, and as techniques 
for extracting sea-surface velocities from the temperature 
data are improved, it is hoped that well resolved surface 
velocity fields will become widely avaible. Analysis of 
these will transform our understanding of dispersal pro- 
cesses in the surface ocean. 
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